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Within the framework of the frozen temperature approximation we develop a strongly-nonlinear 
theory of one-dimensional pattern formation during directional solidification of binary mixture under 
nonequilibrium segregation. In the case of small partition coefficient the full problem is reduced to 
the system of two ordinary differential equations describing the interface motion in terms of its 
velocity and position coordinate. The type of the oscillatory instability bifurcation is studied in 
detail in different limits. For the subcrytical bifurcaton relaxation interface oscillations are analyzed 
analytically and numerically. We show that these oscillation exibit a number of anomalous properies. 
In particular, such oscillations can be weakly- or strongly-dissipative depending on the physical 
| parameters and the amplitude of the strongly-dissipative oscillations is determined not only by the 

form of the corresponding nullcline but also by the behavior of the system for small values of the 
interface velocity. Characteristic parameters of the superlattice occuring in the growing crystal are 
(^i estimated. 
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O ' I- INTRODUCTION 

A great variety of patterns can form during directional solidification of binary mixture when steady state motion of 
the planar solid-liquid interface becomes unstable (for a review see, e.g., J0,||). There are several different mechanisms 
of this instability. One of them, proposed by Mullins and Sekerka pj and called the marginal (morphological) 
, instability, has been investigated and discussed intensively for last years 0)- By this mechanism, under appropriate 
conditions the planar interface loses stability mainly with respect to essentially nonuniform perturbations, which gives 
rise to cellular interface formation. 

Coriell and Sekerka [Q have shown that departure from local equilibrium at the planar solid/liquid interface due to 
high solidification rates can cause oscillatory instability of its steady state motion. Merchant and Davis H analyzed this 
instability in detail taking into account different characteristics of non-equilibrium conditions in a thermodynamically- 
consistent way. Such an oscillatory instability can manifest itself in self-formation of structures parallel to the moving 
interface, which have been observed experimentally for Ag-Cu, Al-Cu, and Al-Fe alloys solidified at high rates |] ||] 
(for a review of materials exhibiting formation of this type structures see |9|). It should be noted that there are also 
other models predicting oscillatory instability of plane solidification front [ 10|-p0[| . 

When solidification proceeds by the normal growth mechanism, non-equilibrium effects at the solid/liquid interface 
become remarkable only at extremely high solidification rates (~cm/sec-m/sec) [Q. In this case within the framework 
of classic approach departure from local equilibrium at the planar interface due to its rapid motion is characterized 
by the interface temperature Ti(v,Ci) depending on the interface velocity v as 

Ti = To + m(v)Ci-- (1.1) 
V 

and the partition coefficient k(v) = C s /Ci is also treated as a certain function of v. Hear To is the melting temperature 
of the pure material being in thermodynamical equilibrium, C s and Cj are the interfacial concentrations of solute 
in the solid and liquid, respectively, r\ is the interface mobility, and m(v) is the nonequilibrium slope of the liquidus 
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regarded as a function of the interface velocity v. Several models for the k(v) dependence at high solidification rates 
have been proposed [^2|-^| and the expression 

= k e + p v 

is much used in theoretical analysis of pattern formation during rapid directional solidification. Here k e is the 
equilibrium partition coefficient usually much less than unity, k e <C 1, and 1/Po is the characteristic velocity being 
typically of the size 5 m/sec The functions m(v) and k(v) are not independent but related by the following 

expression pl| , p5[ 

1 | fc e -fc[l-ln(fc/fe e )] ^ (i 3) 

where m e is the equilibrium value of to. 

In mathematical analysis of directional solidification it is conventionally used the "frozen temperature" approxima- 
tion. In this approximation the latent heat diffusion is neglected and the thermal properties of the liquid and solid are 
taken to be equal, which decouples the temperature and solute fields and enables one to represent the temperature 
distribution T(y) in the frame moving at the pulling velocity V in the direction y in the form 

T(y) = To + Gy, (1.4) 

where G is the temperature gradient at the interface and the origin of the moving coordinate system is chosen so that 
T{y)\ y=0 — Tq. Sivashinsky and Novick-Cohen |^6|,^7j have shown that the difference in magnitude of the temperature 
diffusivity in solid and liquid gives rise to an essentially nonlocal self-interaction of the interface in addition to the 
ordinary renormalization of the temperature gradient G [pi . The effect of latent heat diffusion on rapid solidification 
has been investigated in detail by Huntley and Davis |2q[ , and Karma and Sarkissian p9fl . It has been found that 
latent heat diffusion can suppress the longwave oscillatory instability at high solidification rates. However, because of 
the solute diffusivity D in liquid being much less than the temperature diffusivity Dt ( D <C -Dt) the neighborhood 
of the zero wave number, where the latent heat effect is considerable, is extremely narrow. The latter enables the 
temperature Ti at the interface to be related to its coordinate along the y-axis of the moving frame by the 
expression Q 

m = To+Gm+ ±r <*; «<cot>, (1 . 5) 

Here c p is the specific heat, L is the latent heat of fusion, and (((t)) is the averaged position of the interface at time t. 
When the mean curvature K of the interface is not extremely small so that K ^> 1 but Kh <§C 1, where /it ~ Dt/V 
and h ~ D/V are the characteristic lengths of heat and solute diffusion, respectively, the integral term in expression 



(1.5) can be ignorable. In fact, in this case different interface segments of size 1/K may move independently of one 
another, thus, their total contribution to the integral term should be small, and the latent heat diffusion controls 
solely their averaged position along the y-axis, fixing the value {((t)}. Under these conditions, as follows from 



(1.5), the "frozen temperature" approximation holds also for rapid solidification, at least, in semiquantitative analysis. 



The fact that longitudinal structures forming during rapid solidification are quasi-plane rather than rigorously plane 



is justified not only by the available experimental data but also numerically |30|. 

When crystals grow layer-by-layer, with their surface being practically parallel to one of the singular faces, the solute 
partitioning can become non-equilibrium at slow rates typically used in technological processes, viz., at V ~ 10 -3 -10 -2 
cm/sec The matter is that the crystal growth at the volumetric solidification rate V is attained when atomic steps 
move along the crystal surface at the mean speed V s t C p ~ V/6 CV , where 9 CT is the mean angle of interface misorientation 
from the ideal singular face. For real singular faces forming in crystallization of monocrystals 8 cr ~ 10~ 3 — 10~ 4 [[H], 
thus, V s te P ~ 1 — 100 cm/sec for these values of V. In the layer-by-layer crystallization solute segregation to the 
solid involves two stages: surface trapping of solute atoms by the moving steps and their following incorporation into 
the crystal bulk. For such large values of V^tep the former process becomes essentially non-equilibrium, causing the 
total partition coefficient k{v) to depend on the interface velocity v. As shown by Voronkov |n| in this case the k(v) 
dependence can be represented in terms of 



k e + k s p s (v) 

k{v) = i+m ' (L6) 



where the function /3 s (v) is of the form 
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Ps(v) =01V 



1 . a i ke P 2 
1 + /J 2 v + 



k s (3i 



(1.7) 



the constants j3\, (3 2 are of order 10 2 sec/cm, k e , as before, is the equilibrium partition coefficient, k s is the interface 
partition coefficient, and typically k e <C k s < 1. The given velocity dependence of the partition coefficient for Al 
and P in Si has been observed by Voronkov et al. When the growing interface deviates from the singular face 

the value 9 cr is ordinary about 1° or large and the solute partitioning can become non-equilibrium at substantially 
higher rates of crystallization. A singular crystal surface corresponds to a cusped minimum of the surface energy 
and, thereby, to an extremely large value of the surface tension |33j |. The latter suppresses the marginal instability 
and practically does not affect the oscillatory instability because at such low crystallization rates the effect of latent 
heat diffusion should not be strong enough. In this case the non-equilibrium solute partitioning through the induced 
oscillatory instability can manifest itself in formation of solute bands, where the solute concentration varies only in 
the growth direction perpendicular to the interface, with the solute bands occurring near singular faces only. Such 
pattern formation has been found to form during crystallization of InSb with impurities Se and Te in addition to 
ordinary solute layers caused by temperature fluctuations due to convective instabilities (see p^j35| and references 
therein) . 

There are a number of works concerned with a nonlinear theory of pattern formation in rapid solidification. Braun 
and Davis []36| , [37| and Braun et al. p8|] developed a weakly-nonlinear theory and derived the corresponding Ginzburg- 
Landau equation in the "frozen temperature" approximation. Huntley and Davis p9]] incorporated in this theory 
latent heat diffusion. Novick-Cohen |^7j obtained a Kuramoto-Sivashinsky equation in the limit k — ► 1. Based on 
the theory of oscillations with weak damping |l0| Merchant et al. [[Il| have studied nonlinear behavior of the zero- 
wavenumber oscillatory instability in a small neighborhood of the absolute stability point where dissipation processes 
are not pronounced. The same problem was analyzed numerically by Brattkus and Meiron ]42] ] . Although the found 
phase paths in the i>£-plane corresponding to the interface oscillations deviate from the ellipsoidal shape substantially 
their theory describes actually small amplitude interface oscillations as witnessed by the existence of unclosed phase 
paths in the v£-plane jnj. The appearance of unclosed phase paths in a small neighborhood of the threshold and 
the transition from small amplitude (but may be strongly-nonlinear) oscillations to really relaxation oscillations have 
been analyzed in detail by Baer and Erneux ff3|,[44| having considered as an example the well known "French duck 
solution" of the FitzHugh-Nagumo equation. In the authors of the present paper actually following the boundary 
layer model proposed by Langer et al. for directional solidification have studied strongly nonlinear large amplitude 
oscillations of planar interface and numerically found anomalous behavior of the corresponding limit cycle. It turns 
out that although the corresponding nullcline is of the "N-form" the limit cycle goes near only one of its stable 
branches whereas in the region containing the other stable branch the phase path deviates from it substantially. It 
should be noted that in the phenomenological model for banded structure formation by Carrad et al. the phase 
path has been assumed to follow the nullcline in the conventional way. 

In the present paper we develop a theory of strongly nonlinear relaxation oscillations of quasiplanar solid/liquid 
interface which can occur during directional solidification. The term "quasiplanar" implies that the interface can be 
treated as planar on spatial scales of order h ~ D/V . Since Z?t ^> D and so /it h such an interface does not have 
to be regarded as plane from the heat transfer standpoint, thus, we may assume that the characteristic curvature 
radius R of this interface meets the inequality h <C R <C h^. Besides, as has been mentioned above real longitudinal 
structures forming during directional solidification are not rigorously plane as well as those observed in numerical 
experiments [ [30| . This enables us to make use of the "frozen temperature" approximation, namely, we suppose that 
the interface temperature Tj a nd the interface position £ in the coordinate system moving at the pulling velocity V 



are related by expression (1.4). Dealing with rapid directional solidification this approximation can be justified, at 
least, in semiquantitive analysis especially if we consider such values of the basic physical parameters that do not 
correspond to a neighborhood of the nose of the neutral stability curve, where latent heat diffusion affects solidification 



process most strongly [39 28]. Besides, in this case to a first approximation in the small parameter f/ R the nonlocal 



self-interaction of the interface due to difference in the temperature diffusivities of the solid and liquid vanishes [E6l 



When crystals grow layer-by-layer approximation (1.4) is also justified because nonequilibrium partitioning comes 
into play even at low crystallization rates ]3l| ] at which the latent heat effect is not yet pronounced. 

The second assumption adopted in the present analysis is that the interface remains quasi-planar during oscillations. 
This can be the case when, in particular, the marginal instability is suppressed by large values of surface tension such 
that the absolute stability limit for the marginal instability lies above that of the oscillatory instability. Such conditions 
and the corresponding physical parameters have been analyzed in Q for rapid solidification. For the layer-by-layer 
growth the given assumption is justified to a larger degree because a singular crystal surface corresponds to a cusped 
minimum of the surface energy and, so, to extremely great values of the surface tension |33[ |. The given assumption 
allows us to confine ourselves to one-dimensional equations governing solute distribution in the melt. 

The third basic assumption to be used is the smallness of the solute partition coefficient k(y) for all values of the 
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interface velocity v attained during oscillations. For rapid solidification problem this means that we shall consider 
only the limit of relatively low solidification rates. Concerning with the layer-by-layer gro wth this assumption does 
not impose such exacting requirements. Indeed, in this case, as follows from formula (0|), k(v) <C 1 for all values 
of v if fc e , k s <C 1. The latter inequality is justified, e.g., for such crystals as Si and Ge |31|. Besides, in the present 
analysis we shall ignore the dependence of the liquidus slope m on the interface velocity v because for small values of 
t he partition coefficient k(v) its variations should not lead to appreciable variations of m as results from expression 
(1.3). The assumption of the partition coefficient k(v) being small allows the solute distribution in the melt to be 
treated as quasistationary and enables us to reduce the original problem to a system of two ordinary differential 
equations whose solution can be investigated by the standard methods of nonlinear analysis. As will be seen from 
the results to be obtained below the general properties of interface oscillations are not sensitive to a particular form 
of the k(v) dependence. Keeping the latter in mind, in t he following analysis for the sake of definiteness we shall 
use such a function k(v) that generalizes expressions ( |l.2| ) and ( [l.6|) in a simple way. Besides, for crystals growing 
layer-by-layer and by the normal growth mechanism the velocity d epe ndence of the interface temperature Tj(u) is 
different. Nevertheless, for the same reason we shall use expression (1.1). 

The paper is organized as follows. In Sec. [b] we review the equations specifying the directional solidification prob- 
lem under consideration, briefly discuss linear stability of the steady-state interface motion and recall under what 
conditions the interface may be treated as quasiplanar. In Sec. IB within the framework of the quasistationary ap- 
proximation the full diffusion problem is reduced to the system of two ordinary differential equations. This section 
also includes the bifurcation analysis of interface oscillations in order to compare the results obtained in the quasista- 
tionary approximation with the previous ones. In Sec. IV we investigate the characteristic properties of the interface 
relaxation oscillations that occur providing the bifurcation is subcritical and estimate the corresponding parameters 
of the longitudinal structure forming during solidification. 



II. GOVERNING EQUATIONS. THE COEXISTENCE OF DIFFERENT TYPE INSTABILITIES 

We confine ourselves to the classic model for directional solidification of a binary mixture |l],||,[|j5]] . A more intro- 
ductory and detailed exposition of this process can be found in 0]3Mp6|,[29l and here we only summarize the basic 
equations governing interface motion and our notation. 

The solute concentration C in the melt obeys the equation: 

^=DV 2 C, (2.1) 
where D is the chemical diffusivity in the melt. Far in front of the solid/liquid interface X the concentration C tends 

to Coo 

C — > Coo as y — > co. (2.2) 

The solute diffusion in the solid is neglected. Local conservation of solute at the interface X leads to the following 
boundary condition: 

v n d[l - k{v n )] = -DV n C\ x . (2.3) 

Kere v n is the normal interface velocity, C is the solute concentration at the interface X on the melt side and 
k(v n ) is the partition coefficient. The velocity dependence of partition coefficient k(v n ) due to nonequilibrium solute 
segregation is specified by the expression 

k e + k p*v 

y ' l + f3*v ' y ' 



where (3* is a given constant. This formula is practically just a simple generalization of expressions ( |l.2| ) and ( [Lq ). 

Keeping in mind the aforesaid in Sec. | we confine our consideration to the "frozen temperature" approximation in 
which the temperature 7$ at the interface X and its y-coordinate C, in the frame moving at the pulling velocity V in 
the y-direction are related by the expression 

T = T Q + G(. (2.5) 

Bere To is a reference temperature chosen to be equal to the solidification temperature of pure melt (C = 0) and G 
is the imposed temperature gradient. 
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Following we suppose that the deviation of the interface temperature Ti from the local equilibrium value T cq 

causes the crystal growth at the velocity: 

V n = r/( T cq - Ti), (2.6) 

where rj is the interface mobility assumed to be isotropic and constant. The local equilibrium temperature T eq , solute 
concentration d at the interface X and the interface curvature H are related by the Gibbs-Thomson condition: 

T e(l = T (l-dH)-mC i , (2.7) 

where m > is the slope of the liquidus line in the binary phase diagram, and d is the capillary length which is 
proportional to the surface tension and accounts for anisotropic effects. The value of m may be treated as constant 
for such v n tha t k( v n ) <C 1. In paper || the interface mobility rj has been represented actually as ry = Vo/m. 



Expressions (2.1)-( [2.7| ) form the comple te m o del for the directional solidification process under consideration. 



The steady-state solution of equations (2.1)- fl2.7|) describes the motion of the planar interface X at the constant 
velocity v st — V, the uniform distribution of impurity in the solid C|*j = Coo, and the stationary solute distribution 
in the liquid , y > C 

C st =C ao (l+ (j±-r - 1) exp[-(y - ■ (2.8) 

Under certain conditions the steady-state motion of the interface loses stability and two different type instabilities, 
oscillatory and marginal, can occur ||f|]. The coexistence of these instabilities has been analyzed in general in 
P,p| j29|j39[ | . However for the purpose of the present paper it will be the best to consider briefly this problem in the 
limit of small partition coefficient (k -C 1). 

As follows from the results to be obtained below the amplification rate a of a small perturbation £(x,t) ~ 
cxp[j}(aVt + ipx)] of the planar interface (where p is the dimensionless wave number) can be considered to be 
substantially less than one, |aj<C 1, as long as k(v) <C 1. In this limit the standard linear stability analysis leads to 
the following eigenvalue equation: 

s J (l + /i) + u9(p) + »(p)=0 1 (2.9) 

where the functions 

Z(p) = ^+ t ik-Vk' v + (2 + 7 + / i)p 2 , 
»(P) = ^ + ljj+jk-l]p 2 + 7P 4 , 

and the parameters 

VtoCoo T dVk Vk 



M 



GDk DmCoo r\m,C a 



It should be noted that Merchant and Davis |^] used practically the same parameters where, however, the partition 
coefficient takes the equilibrium value k e . 



According to equation (2^) the planar interface becomes unstable (i.e. Recr > 0) with respect to the given 
perturbation when SR(p) < or 3(p) < 0. The former inequality is associated with the marginal instability, which 
is caused by the constitutional undercooling and has been thoughtfully analyzed and discussed (for a review see 
p|,pPj4^| ) . The corresponding maximum of the amplification rate cr max is attained at p = p max 7^ and, thereby, the 
neutral curve I, separating stable and unstable (from the marginal instability standpoint) regions in the CooV'-plane, 
depends on the capillary length d (Fig. Q). This instability region is specified by the inequality 5i(pmax) < 0, leading 
to the expression 

DG kT dV 



mCoo > + — ^— + 2y/T Gdkj. (2.10) 

The latter inequality is bound up with the oscillatory instability which is described in |4|,|| and occurs through the 
dependence of the partition coefficient k(v) on the velocity v. Since the function is increasing the corresponding 
maximum of amplification rate a is attained at p — 0, thus in the "frozen temperature" approximation the oscillatory 
instability region in the CooV - plane is determined by the inequality 
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and is independent of the capillary length d. The latent heat diffusion can suppress the oscillatory instability for small 
wave numbers p9| , p9[ |. However due to the heat diffusivity being much larger than the solute one the wave number 
interval where the latent heat essentially affects perturbation dynamics is narrow enough. Therefore, the condition 
( 2. 11| ) obtai ned i n the "frozen temperature" model holds true if the latent heat fusion is not extremely great. The 



equality in (2.11) specifies the n eutra l curve II of the oscillatory instability (Fig. 

In particular, it follows from ( [2.11 ) that for arbitrary values of the pulling velocity V and the imposed temperature 



gradient G all disturbances corresponding to the oscillatory instabiblity decay when Coo < C* , where 

k 2 

(j* — 3 (2 12) 

mq{k s - fc e )/3* ' 

In other words, the inequality Coo < C* is the absolute stability criterion for the oscillatory instability in the limit 
of small partitition coefficients. It should be noted that this criterion matches the results obtained by Merchant and 
Davis H in the limit V — > 0. For > C* the oscillatory instability can occur only in the region > C as (/3) 
specified by the function C as (/3) of the dimensionless pulling velocity (3 = j3*V 

C as ({3) =C*(1 + ^) 2 , (2.13) 
r 

where the parameter r — k e /k s . In the Coo^-plane the curve Coo = C as ([3) is the boundary of the absolute stability 
region (line 1 in Fig.|). In this fi gure line 2 represents the neutral curve of the oscillatory instability for a fixed value 
of G measured in units of G* = k e j (Drj/3* 2 ). In particular, in these terms the neutral curve can be represented as 



Coo — C as (/3) 



G 1 + /3 



(2.14) 



G* l3 2 (l + /3/r)_ 

The neutral curve in contrast to the absolute stability curve attains its minimum at V — V cyi t > 0. For example, as 



follows from (|2T4j) when r < 1 and r 2 « G/G* « 1/r the value V cxt belongs to the region r 2 / (3* < V <C and 
so 

V«^Q§). (2.15, 

As the capillary length increases the neutral curve I for a fixed value of G goes upwards in the Coo^-plane and, 
thus, there can be three characteristic relative positions of curves I and II shown in Fig. |l[ When the capillary length 
is small enough (d <C di, where d\ = di(D, G, r/, To) is such a value of the capillary length that curve I crosses curve 
II approximately at its minimum point) the region of the oscillatory instability is totally located inside that of the 
marginal instability. In this case the oscillatory instability of the planar interface seems not to be able to manifest itself 
in any way because of rapid development of spatially nonuniform interface perturbations (Fig. 0(a)). When d ~ di 
there is a certain small region in the Coo^-plane that corresponds solely to the oscillatory instability (Fig. 0(b)). In 
this case the oscillatory instability can occur. However, as it has been shown numerically pC| j, the interface becomes 
substantially nonplanar during the instability development. In this case the banded structures seem to occur [||. For 
d d\ the large value of the capillary length, i.e. of the surface tension, suppresses the marginal instability and 
thereby the interface tends to be planar. It should be noted that in this case (Fig. ||(c)) the developed form of the 
solidification front will be quasiplanar and solute band formation is expected [gj. 



III. QUASISTATIONARY APPROXIMATION 

In this section we obtain nonlinear evolution equations for the planar interfacial position y — (,{t). In general, these 
equations must contain memory effects; that is a displacement of the interface X from its steady-state position causes 
a perturbation of the impurity distribution which, in turn, affects the motion of the interface at later times. However, 
in the limit of small partition coefficient it turns out that the interface velocity varies in time so slowly that it remains 
practically constant during the time needed for relaxation of the solute distribution in the melt. In other words, it 
means that all temporal scales of interface dynamics are much larger than the characteristic time (t v ~ D/V 2 ) it 
takes for the steady state solute distribution to form when Ci and v are fixed. 



G 



In the frame of reference moving in the y-direction at the interface velocity v the diffusion equation has the form: 



C — DCyy ~f~ vCy. 



(3.1) 



Following [|4G|j4!j|3(J we integrate (3J) over y from £ to oo then taking into account the boundary condition ( p.3| ) we 
obtain 

^ poo 

(3.2) 



- / dy[C{y,t) - Coo] = «[Coo - k(v)d]. 



In the quasistationary approximation we can regard the left hand side of equation (^J) as a small pert urb ation and 
neglect it at lower order of the perturbation technique. In other words, setting the transient term in (3.1) equal to 
zero and then solving the obtained equation we find the lower approximation of C(y,t): for y > C 



an.h = (\ [Gm-r^rxi.-j-^fo-o 



Substituting expression (3.3) into equation (3.2) we get 



D- 



dt 



Ci Co 



^[Coo - k{v)Ci 



Equation (3.4) along with the equation 



and the relationship 



— = v — V 
dt 



v = -r](G( + mCi) 



resulting from ( |2.5| )-(2.7) form a complete description of the quasiplanar interface motion. 
The steady state solution of the given system of equations is of the form: 



Co 



k(vy 



c st = - 



G 



V mC, 



k 



It is convenient for the most of the following analysis to use the dimensionless variables, a and u: 



c-c st 

a = —pi-> u 

In terms of this variables, the equations of motion become: 

du (l + u) 3 k[u] 



dt 

da 

t—t = 

dt 

Here we have also introduced the quantities: 

D 



(l + o) k[0] 
— (tuj) 2 u. 



v-V 
V 



a [u}) 



i 



< k[u] = k(V [l + u}), 



V 2 k ,y ' ' kM(l + /.*)' 
where u> is the frequency of the interface oscillations at the threshold, and the function 



1 



a it 



(1 + /i) \k[u] 



+ flU — 1 ) + (to) 



k[u) (1 + u) 2 



(3.3) 

(3.4) 

(3.5) 
(3.6) 

(3.7) 

(3.8) 

(3.9) 
(3.10) 

(3.11) 
(3.12) 



specifies the nullcline of equation (J^ 

In this way the full diffusion problem of interface dynamics is reduced to the system of ordinary differential equations, 
which describes the nonlinear oscillations of the quasiplanar interface Ea] . 
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In the present section on the basis of evolution equations ( |3.9| ), ( [3.1C| ) we also study the bifurcation mode of the 
oscillation onset, depending on relation between physical parameters such as V,G,rj etc. In addition, to validate the 
quasistationary approximation, we compare the form of the interface auto-oscillations obtained in the quasistationary 
approximation for the supercritical bifurcation and those obtained in the one-dimensional full diffusion problem 
described in Sec. [n| 

In order to analyze the bifurcation mode we can make use of the Bogolyubov-Krylov-Mitropolskii technique and 
the theory of Hopf bifurcation |pl[ , p2| . In this way near the threshold the interface motion is described by the 
quasiharmonic time dependence of the variable a = A(t) cos(uj(t)t), where Lo(t) ~ u and A(t) are certain functions 
of the time t being practically constant on the temporal scale 1/ui. Following the standard procedure |^l],[52| we find 
that the amplitude A(t) of the quasiharmonic oscillations obeys the equations 

d A 1 

r— = -A(e + aA 2 ), (3.13) 



where the constants e and a are specified by the formulae 

(l + /i)fc[0] (.'■' M 



a o[0] = n , Lrm (W] + i - k'[0] ) , (3.14) 



8(rw) 2 



C[0] + 2oJ[[0](3- 



k[0] 



(3.15) 



and the primes on the symbols denote the derivatives of the corresponding functions with respect to the variable u 
taken at the point u = 0. In these terms the instability condition of steady state interface motion takes the form 
Oq[0] < 0. In other words, for the oscillatory instability to occur the n ullclin e a = ao{u] should be decreasing in the 
vicinity of the steady state point (u = 0, a = 0). Taking into account ( 3.1 2| ) we can rewrite the latter inequality in 
the form /i + l/(fc[0]Af) < /c'[0]/fc[0] which as it must is exactly the dimensionless form of inequality (2.11). 

In the following analysis it will be convenient to measure the parameters /i and l/(fc[0]M) also in units of fc'[0]/fc[0]. 
In other words, let us introduce new parameters (fi g , <fi^ L > defined by the formulae 



k[0]M 



9 m ■ 



M = <t>. 



M fc[0] ' 



(3.16) 



In particular, at the points of the neutral curv e <fi a + (fi^ = 1 an d at the absolute stability boundary, /i as = fc'[0]/fc[0], 
we get <fi g = and (fi^ — 1. Then substituting ( 3.12 ) into ( |3.15| ) and calculating the obtained result at the thre shold, 
a[j[0] = 0, we find an expression for the parameter a which for the k(v) dependence specified by formula (2.4) takes 
the form 



a = a [u] (4> - (fi b [u\) . 



(3.17) 



Here the parameter (fi = <fi g G [0, 1] characterizes deviation of the system from the absolute stability boundary, ao[ 
is a certain positively definite function of the variable u, and the function <fib[u] is given by the expression 



J r + fil 1 + /3 r + /3 

(l-r)/3 ,„ 2(3 
3 + - — (4 - 



{l + (3){r + (3Y 1 + p r + /3 



(3.18) 



According to (3.18J ) bi furcation of the oscillatory instability is supercritical if <p < <fib[P] and subcritical for tfi > <fib[0\- 
As follows from (|3.18| ), when k e <C k s , i.e. r <C 1 there are three limits characterizing different behavior of the 
function <fib[0\- In the limit j3 -C r (i.e. V <C r/(3*) the function <fib[/3] w 2/3/r which actually exactly matching 
the results obtained by Huntley and Davis |3j| for small pulling velocities in the frozen temperature approximation. 
When r < (3 < 1 (i.e. r/j3* < V < l/f3*) the value <f> b [(3] w 5/6 and for ~ 1 (V ~ 1/0*) we get that (fi b [0] -> 1 
as (3 — > 1 and cj) b [0] > 1 when (3 > 1. In the latter case, > 1, bifurcation is always supercritical. In the CooV- 
plane the curve corresponding to the transition point from supercritical to subcritical bifurcation is specified by the 
formula C tr = C ajS (0)/(l — <fib[/3}) (for j3 — f3*V < 1) and is shown by the curve 3 in Fig. ||. For a fixed value of the 
temperature gradient G the neutral curve in the C^y-plane crosses the curve Coo = Ctt(0) at a certain point (V c , C c ) 
and for V < V c as well as V > V c the oscillatory instability onset at the threshold is characterized by subcritical 
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and supercritical bifurcation, respectively. In particular, when r 2 <C G/G* « 1/r this point belongs to the region 
r//3* < 7 < 1//3* and in this case the values V c , V cxt are related by the expression V c = (2/5) 1 ' 3 V e ^t- 

In the region V > V c where bifurcation is supercritical and the coefficient a is negative the amplitude of the 
steady-state quasiharmonic oscillations Aq near the threshold is 



oc 



f 



m 

fc[0] y k[0}M 



-,1/2 



/') 



(3.19) 



For V < V c , the coefficient a is positive and in order to find the amplitude of the interface oscillations a more 
complicated analysis is required (this case is considered in Sec. IV). 

It should be noted that the results of the present section concerning the bifurcation behavior are in qualitative 
agreement with those obtained in the previous papers for different physical conditions |l^ , |3^ -^8| . 

To examine the feasibility of the quasistationary approximation in the supercritical mode we compare the lim it 
cycles of the interface autooscillations as well as t he ti me course a(t) and u(t) obtained by solving equations (|3.S|). 
(3.10) and the full one-dimensional problem ( p.lj )— (2.7). For this purpose we have solved numerically the equations 
mentioned above using the following values of the parameters: fco = 0.0125, fcoo = 0.8, (3*V = 0.125, kM = 10, and 
/i = 0.5. In this case the ratio fc'[0]/fc[0] ~ 0.75 and, thus, in the C^F-plane the equilibrium point of the system is 
practically located near the boundary of the oscillatory instability region because of [k'/k — (k/M + /i)]/ [k 1 / k) ~ 0.2. 
The obtai ned re sults are plotted in figures |]-[| where the dashed and solid lines corresponds to the solutions of 
equations (|3.9|), fl3.10t ) and the full diffusion model, respectively. 

As seen from Fig. |3j the quasistationary approximation gives an adequate description of the interface oscillations at 
least in the supercritical mode. Fig. ^ also demonstrates that even in the vicinity of the instability region boundary 
the nonlinear effects play a significant role in the interface oscillations and the limit cycle substantially deviates from 
the elliptic shape and the time dependence u(t), a(t) are practically of the spikewise form. Nevertheless, the period 
T of these interface autooscillations can be estimated by the expression T — 2-kuj^ 1 — 2ttt[{\ + ^)fcM] 1//2 which is 
rigorously justified for quasiharmonic oscillations only. Fig. [5] shows the resulting impurity distribution in the solid 
C so \(y). According to Fig. || and as would be expected the spatial period of the growing superlattice is of order 



H ~ TV ~ 2tt 



(1 + fi) 



Vk 2 G 



1/2 



(3.20) 



IV. RELAXATION OSCILLATIONS OF INTERFACE 



As shown in Sec. Ill the subcritical bifurcation corresponds to small values of the pulling velocity, namely, V < 
((3 < 1). Besides, the asymptotic behavior of the curve separating in the CooV-plane the regions of the subcritical 
and supercritical bifurcation is different for (3 <§; r and r <§; j3 <§; 1 (when r = k e /k s <C 1). So for (3 < 1 the interface 
oscillations may be expected to be relaxation, with their main features depending substantially on the ratio of (3 
and r. In agreement with the results to be obtained below there are two cases of strongly-nonlinear dynamics of 
the interface motion: weakly- and strongly-dissipative oscillations. Weakly-dissipative oscillations (however of large 
amplitude with respect to u) can occur when 

(3 « r and - « (1 - 9 « ( -) ' . (4.1) 



r 



whereas the conditions 



1/3 

(3 < r and I - ) < (1 - (j) g (4.2a) 



r < (3 < 1 (for r = k e /k s < 1) (4.2b) 

will result in strongly-dissipative oscillations. It should be noted that in the strict sense only strongly-dissipative 
oscillations are relaxation, although in the two cases the velocity amplitude u can attain large values much greater 
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than unity. In this section we analyze strongly-nonlinear dynamics of the interface oscillations in the given cases 
individually. 

Weakly-dissipative oscillations of crystal interface that can occur during directional solidification has been analyzed 
in detail by Merchant et al. Ell and Brattkus and Meiron 0] within the framework of the classic model for rapid 
solidification similar to one used in the present analysis. So we only briefly consider the first limit (4.1) in order to 
complete the analysis of strongly-nonlinear dynamics of the given model. 

In this case the ratio 



k'[0] (1 - r)f3 ^ 

k[Q] (l+j9)(r + j9) ~ r 



is a small parameter, and the value 



(™) 2 ~4>g~- 



According to Sec. [[TJ, for (3 -C r the bifurcation is subcritical when tj> g > 2(3 /r and at the threshold (j> g + (j)^ = 1. 
So in order to analyze relaxation oscillations of the interface in the given limit we may confine ourselves to the case 
1 — </>n ~ 4>g an d <t> a ^ > (3 jr. We note that the two relations have actually led us to the le ft-h and side of the latter 



conditions of limit (4.1). Its right-hand side enables us to treat the ratio ao[u]/a in equation (3J3) as a small parameter 
for the interface oscillations that will occur under these conditions. The matter is that the characteristic value of 
the amplitude a of the variable a turns out to be about a ~ toj, so a ~ (cjigP/r) 1 / 2 and the maximum u max of the 
interface velocity u attained during oscillations can be estimated as u max ~ w m i n , where Umin is the u-coordinate of 
the point ({t m j n , Q m j n ) at which the curve ao[u] attains its extremum (minimum) in the region u > (Fig. |^). In this 



case expression (3.12) may be rewritten in the form 



M J 
r 



-(i-^) + 7« + ^p^5 



(4.3) 



whence for cf> g *~ 1 



b„ we get that 



v v 

Umin ~ ^(1 _ 4>n) ~ > 1, 



2{3 



/3 1 



(4.4a) 



Imin ~ -7(1 - </V) Z 



(4.4b) 



so, by virtue of (f4.l| ), we obtain that |ao[w]| ^ |a m i n | a. 

At lower order in the small parameter ao[u]/a the system of equations ( p.9[ ), ( |3 .10 ) is conservative with the "energy" 
(the first integral) 



H(u, a) 



(™) 2 



(1 



1 2 

—a 
2 



(4.5) 



When deriving expression (4.5) we have also taken into account that for the given values of the parameters k[u] « k[0] 
and a <C 1. The remaining term in equation (3.£) proportional to a[u] causes time variation of the "energy", namely, 



^ I \2 r i 

r— = -(Tui) a [u\u. 



(4.6) 



In the case under consideration the system dynamics may be described in terms of fast motion along a phase path, 
specified by the equation TL(u,a) = h for a fixed value of h, and slow time variations in the "energy" H.(u,a). The 
geometry of the phase paths in the ua-plane is shown in Fig. |(]. As seen in Fig. ^ there are two types of the phase 
paths, closed and unclosed, which correspond to h < i(™) 2 and h > ^(tui) 2 , respectively, and under the adopted 
assumptions solely the closed paths are meaningful. 

Steady-state oscillations matches such a value h s t for which the right-hand side of equation (|4.6| ) averaged over 
a single period of the system motion is equal to zero. Since the right-hand side of equation ([4.6j) depends on the 
variable u only, at lower approximation we may average it assuming the system to move strictly along the phase for 
a certain fixed value of h. In this way from (4.6) we get 
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dh 2t 



<(h) 



dt T(h) 



da + {h,u) 
du a ° ^ 



(4.7) 



where T{h) is the period of the system motion along the phase path corresponding to 7i = ft, it m in (ft) and it max (ft) 
are the minimum and maximum of the variable u attained during the motion along this phase path, and a + (h,u) is 
the positive solution of the equation Ti{u, a) = ft, namely, 



Umin(M 



(2h) 



1/2 



(2ft)V2 



^max(^) 



(2ft) 



1/2 



tui - (2k) 1 / 2 



(4.8) 



and 



a+(ft, u) = (™) 



2/i 



1 + it 



1/2 



(4.9) 



Substituting (4.8) and ( f4.9| ) into (4.7) and integrating over it we obtain 



dh 

T —— — TT(TUJ) 



/3t 2h 



dt " v '" / rT(ft) (rw) 2 
2/i 



x (1 - 0„)Ai 



(tw) S 



<?A 2 

r 



2h 



(4.10) 



where the functions 



Ai[x] = 



vT^¥(l + y/T=x) ' 



(4.11) 



A 2 [x] = 



2x(\ + 2Vl T7 x) 



(l-a;) 3 / 2 (l + \/T 3 £) 2 



(4.12) 



In agreement with the results obtained in Sec. [II from (4.10) we can see that the oscillatory instability occurs when 



< 



1 and is subcritical for </> M + 2(3 /r < 1. As follows from ( 4.10| ) the dependence of the "energy" /i st of steady 



state oscillations on the physical parameters is specified by the expression 



4> g = (1 - </> M )Ai 



2h s 



A 2 



2ft s 



(4.13) 



In particular, for (3/r < 1 - < (/3/r) 3/2 and 0„ < [(1 - Al )(r//3) 1/3 ] 3/2 the value 



(™) 2 



/3 



(4.14) 



Besides, according to (B^q) the maximum a max and the minimum a m j n attained by the variable a during oscillations 
~ (2ft, s TJ f/2 . Whence for the given values of <\> g and M we obtain 



^max ~ u max 



^max ^ "75" (1 0/0; ^ n 



(4.15) 



The period T of these oscill ation s is practically determined by the time r s it takes for the system to pass through 
the region u < 0. Equation (3.10) enables us to estimate this time as t~ ss Ta max / (tlo) 2 whence we get 



T ~ — . 



(4.16) 



In the case under consideration the characteristic relative position of the nullcline ao[u] and the limit cycle of the 
oscillations are demonstrated in Fig. 0. 
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In limits (4.2a) and (4.2b) the interface oscillations are strongly-dissipative. The corresponding phase paths as a 
whole cannot be described by level curves of any function similar to the energy H.(u,a). In particular, in the region 
u < 1 these phase paths may be roughly related to the unclosed level lines of the energy Tt(u, a). In order to analyze 
the interface oscillation s in t his case we may make use of the classic theory of relaxation oscillations ]53f| . First, we 
consider in detail limit (4.2a). Under these conditions, as it can be shown directly from the system of equations (|3.S|), 
( |3.10 ), the phase path going through a point that is not located in the region \u + 1| -C 1 or in a small neighborhood 
of the nullclinc a = ao[u] forms practically a horizontal line, at least, in the vicinity of this point. In other words, 
at such a point the value du/dt is large in comparison with da/dt and the quantity tlo may be treated as a small 
parameter. Therefore, the interface dynamics should include fast motion governed by equation (3.9) for fixed a unti l 
the phase path reaches a small neighborhood of the nullclinc a = ao[u], where the left-hand side of equation (3.9) 
tends to zero, or a small neigh borh ood of the boundary u = — 1, which formally can be also regarded as the other 
branch of nullcline of equation ( |3.9| ). In the two latter regions the interface motion is slow. 

The nullcline a = ao[u] (shown by curve 1 in Fig. ^) looks like an "N" and the stationary point (u = 0,a = 0) is 
unstable when it belongs to the decreasing segment of the nullclinc. So according to the classic theory of relaxation 
oscillations |^3| one could expect that th e li mi t cyc le will be of the form presented by the dashed line in Fig. [?|. 
Indeed, based on the system of equations (3.9), (3.1C) we can show that the increasing segments I, III of the nullcline 
a = ao[u] are formally attractive for the phase path whereas the increasing one is (segment II) is repulsive. The 
main characteristics of such a limit cycle are directly determined by the form of the nullcline. In particular, the 
minimum and maximum attained by the variable a during oscillations are approximately equal to a m ; n and a max , i-e- 
the a-coordinates of the extremum point s (u m i D ,tt m i n ), (u max , a max ) of the nullcline. In the case under consideration 
the curve ao [u] is specified by expression ( |4.3| ) and the last term on its right-hand side is ignorable for ii» 1. Whence 
we find that the quantities a min , u m i n can be estimated by expressions (4.15) and for a max , u max we obtain 



< 1 < a max < 



(1 - 



(4.17) 



However, the right-hand side of equation (3.9) possesses a certain peculiarity in the region \u + 1| <C 1, because it 
tends to zero as u — * — 1. Therefore, if the phase path goes into this region the time scale hierarchy will be changed 
and, thus, the phase path will not be able to follow the segment of the nullcline a = ao[u] going through this region, 
and the maximum a max attained by t he va riable a during oscillations will be not determined by the extremum point 
(w max , a max ). This is the case in limit ( 4.2a ) because for such values of the physical parameters the nullcline a = ao[u] 
reaches a small neighborhood of the line u = —1 in the region a > and the phase path has to go into the region 
\u + 1| <C 1 returning from large values of u. Therefore, under the given conditions the limit cycle of the interface 
oscillations will involve the conventional segments going along the horizontal lines a ~ a max , a — a m i n , along segment 
III of the nullcline a — ao[u] from a max to a m j n in addition to a certain anomalous segment located in the region 
\u + 1| <C 1 which joins the two horizontal segments. 

In order to complete the limit cycle construction we need to analyze the system motion in the region \u + 1| <C 1 
by solving directly the system of equations ( p.9[ ), (3.1C). Dividing equation ( |3.9| ) by equation (3.10), setting u = — 1 
except for the terms containing the cofactor (1 + u), and taking into account that in the given case a<Cl we obtain 
the following equation describing the phase path in the region |1 + u\ <C 1 



du 
da 



(1 + uf 



1 



(1 



1 



(4.18) 



subject to the formal "initial" condition 



(4.19) 



Condition (4.19) reflects the fact that the limit cycle segment located in the region |1 + u\ <C 1 originates from the 
phase path going into this region practically along the horizontal line a = a m i n . To second order in a the solution of 
equation (4.18) meeting ( 4.19| ) is of the form 



1 r 

JZ j ~ „ , (a — 5.min)( — a — Omin)- 



(4.20) 



As it must u — > o o as a — > a m i n + 0, the second point where the value u tends to infinity is a = — a m i n - Besides, 
according to ( 4.20 ), values of u such that (1 + u) ~ 1 correspond to values of a belonging to small neighborhoods of the 
points — a m j n and a m ; n , i.e. |a + d mul \ <C |a m i n | or |a — a m i n | <C |a m in| because in the given case the ratio ra^ in / (/30 9 ) 
is a large parameter. Whence it follows that the desired value of a max can be estimated as a max ~ — a m j n 



min/ 

i.e. 
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(4.21) 



It should be noted that a max 3> a max , so the limit cycle goes over the externum point (u m ax, a max ) at a remarkable 
distance. In addition, formula ( 4.20| ) shows that the minimum u m i n attained by the variable u during oscillations is 
about u min w -1 + [(/3(/> s )/(ra^ in )J i/2 , thus 



-1 



-,1/2 



r(l - M ) 4 



(4.22) 



The maximum Umax attained during oscillations is determined by the intersection point of the horizontal line a = a,, 
and the nullclinc a = ao [u] , whence we get 



(l + V2)(l-^)r 
2(3 



(4.23) 



Let us now estimate the characteristic time scales of the interface oscillations. According to the classic theory of 
relaxation oscillations |53f| the fast motion along the horizontal lines a — a m i n and a — a max governed by equation (3.E) 
is characterized by the time scale Tf ~ T/a max . The time during which the system moves slowly inside the region |1 + 
u\ <C 1 practically in the a-direction from the point a m i n to a max is about t~ ~ 2a max r/ (tuj) 2 ~ r(l — M ) 2 r/ (/3<p g ) 3> 
r/. The motion along segment III of the nullcline a = ao[u] is actually governed by equation (3.10) where we may 
roughly set u ~ u max . In this way we find that for this stage of the interface motion the characteristic time can be 

r. It should be pointed out that in spite of r+ <C Tf ~ T/a max the 



estimated as r+ ~ r s /u max ~ r(l - 4>ii)/4>g 
motion along the nullcline ao [u] in the region u 
lines a = a max ora = a min 
given case the period T of the system motion is determined, as before, by the value r, 



is slow in comparison with the fast motion along the horizontal 



because in this region the time scale of the fast motion is Tf / u^ ax rather than 



thus 



T f- 



r (1 - <j> T f 

1 ~ T 

& <\>S 



In the 



(4.24) 



(Sec. Imp 

nullcline cio[u 



Let us now consider the interface oscillations in limit ( [4.2b ). In addition we shall assume that 4>^ -C </> M & = 1/6 
because for <j)^ < (p^ the oscillation onset is subcritical. In this case the k[u] -dependence and the 

may be approximately specified by the expressions 



k[u] » fc[0](l+u) 



(4.25) 



and 



a [u\ 



1 



1 



(4.26) 



For (j) g ~ 1 the system of equations (3J5), ( 3.10 ) does not contain actually a small parameter. Nevertheless, the 
minimum a m ; n attained by the nullcline ao [u] in the region u > is located in the immediate vicinity of the boundary 
a = — 1, i.e. a m ; n ss — 1 and matches the value £t m i n ~ l/(</> A1 ) 1 / 2 3> 1. Besides, the increasing segment of the nullcline 
ao[u] (for u > 0) corresponds to large values of the variable u. Therefore, in the given case the variable u attains 
large values during interface oscillations and the limit cycle comes close to the line a = — 1. Since the right-hand 

) 4 /(l + a) this behavior of the limit cycle causes the interface 



side of equation (3.9) depends on u and a as (1 



oscillations to be relaxation and the limit cycle contains segments of the fast motion that approximately are parallel 
to the it-axis, a segment following the increasing part of the nullcline a = ao[u] in the region !i» 1, and a segment 
located in the region \u + 1| < 1. The latter one, however, also approximately follows the corresponding increasing 
part of the nullcline a = ao[u] as results from numerical analysis although the bifurcation is subcritical for the given 
physical parameters. This characteristics explains the fact that the region of the parameter <p g where the interface 
motion is linearly-stable and the interface oscillations with a finite amplitude can occur is sufficiently narrow. In 
particular, as formally </> M — > this region is determined by the inequality 1 < <f> g < (j> gc pa 1.05. 

Inside the instability region far from the neutral curve, i.e. for <j) g <C 1, the limit cycle again will reach the region 
\u + 1| -C 1, thus the corresponding segment of the limit cy cle w ill not follows the increasing part of the nullcline 
a = ao[u] and take a form similar to that matching limit (4.2a). This conclusion has been obtained by solving 
numerically the system of equation (3.9), (3.10) in the region |u+ 1| <C 1 for <f) g <C 1 where it can be rewritten in the 
parameterless form 
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du p a 



where u v 



(u- 



-l)/^ /2 , 



dt 
da p 



a p 



Up 



(1 



(4.27a) 
(4.27b) 



= (a+ l)0g 7 ^, and the time t is measured in units of r/ifig'*. In the given region the phase 
path substantially deviates from the nullcline and goes over its extremum point (u ma xj Omax) at a certain distance. 
So also in the case under consideration the limit cycle of the interface oscillations and one formally predicted by the 
classic theory of relaxation oscillations [|53| differ significantly in form. The period of such oscillations can be roughly 
estimated as 



T 



6f 



(4.28) 



Concluding the present section we note that, first, in all the cases considered hear the interface oscillations may 
be regarded as spikewise. T he m at ter is that the interface dynamics described in terms of the variables a and u, i.e. 
by the system of equations ( |3.9|) , ( |3.10|) is characteriz ed b y self-acceleration in the region u -C 1, which is reflected 
in the existence of the cofactor (1 + u) 3 in equation (3.9). So the interface dynamics is characterized, at least, by 
two time scales, t~,t + , corresponding to the motion in the regions \u\ < 1 and u» 1, respectively, with t~ r + . 
Therefore, the time course of these oscillations contains spikes of duration r + formed during motion in the region 
ii> 1. Second, from the standpoint of the crystal growth theory it is important to analyze the main characteristics 
of impurity distribution in the growing crystal. So we also estimate the spatial period H of the impurity distribution 
C so \{y) occurring during strongly-dissipative interface oscillations. By definition 



T„ 



V / dt{l + u). 



(4.29) 



For example, in limit (4.2a) the main contribution to the value of integral ( 4.29Q is d ue to the system motion along 
segment III of the nullcline. Then substituting the corresponding estimates into (4.29) we get H ~ l^Mmax^, so 



H 



Vk(V) Pcbg 



In case ( 4.2b ) for <p g < 1 and <C 1 the spatial period is about 



H 



D 



Vk{V) 



(4.30) 



(4.31) 



The comparison of ( |4.30|) and ( 4.31 ) shows that expression ([4.31 ) can be treated as a rough estimate of the spatial 
period of the impurity distribution occurring during strongly-dissipative interface oscillations. Keeping in mind the 
layer-by-layer crys tal g rowth we set D ~ 5 • 10~ 5 cm 2 /sec, V ~ 10~ 2 cm/sec and k(V) ~ 0.1 for such values 
of V. Then from (4.31) we get H ~ 500 fim, which coincides in order with the characteristic length of impurity 
microinhomogeneities observed in real crystals . 

In addition, we consider some characteristic features of the C so \(y) dependence. Expressions (3.6), fl3.8| ) lead us to 
the following relationship 



a 



sol 



= 1 



k[u] . 
-^r[l + (1 + fi)a - fiu\. 



(4.32) 



During the system motion practically alon g segment III of the nullcline do [u] the concentration C so \ is actually equal 
to Coo beca use t he right hand side of ( 4.32j ) differs from one by the value (tlu) 2 u/ (1 + u) 2 -C 1 for u ^ 1. As seen from 
expression ( [1.32] ) in cases where the impurity concentration in the solid near the interface attains its minimum 

C™J n and maximum C™f x when in the ua-plane the system reaches the points with the coordinates 
and a ~ a m ax, u ~ w max , respectively. So the minimum and the maximum of the impurity distribution C so i(y) can be 
estimated as 
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sol °c 



k[0] 

k " ma 



■(1 + a m in), 



In order to analyze the feasibility of t he q u asista tionary approximation for the interface rel axat i on o scillations we 
have also solved the system of equations (|3.9| )-( 3.1C ) as well as the full one-dimensional model (2d)- (2/7) numerically. 
We used the following values of the parameters: fco = 0.0125, fcoo = 0.8, (3*V = 0.125, Mk = 10, and /i = 0.2, which 
correspond to the subcritical bifurcation and are not too far from the critical point. The obtained results are presented 
in Fig. §-Fig. [l(]. 

It should be expected that the quasistationary approximation may be violated when the interface velocity becomes 
small enough: (1 + u) <C 1. Nevertheless, as seen from these figures, the shapes of the limit cycle constructed by 
different methods qualitatively coincides with each other. Therefore, first, based on the quasistationary approximation 
one can describe the main characteristics of oscillatory zoning. Second, the estimates for the amplitude and spacing 
of the growing superlattice can be used as a first approximation for the real processes. 



ACKNOWLEDGMENTS 



The authors thank Professor V.V. Voronkov for interest and useful discussions. The research described in the 
present paper was made possible in part by Grants No. U1I000, U1I200 from the International Science Foundation. 



[1] J. S. Langer, Rev. Mod. Phys. 52, 1 (1980). 

[2] J. D. Weeks and W. van Saarloos, J. Cryst. Growth 112, 244 (1991). 
[3] W. W. Mullins and R. F. Sekerka, J. Appl. Phys. 35, 444 (1964). 
[4] S. R. Coriell and R. F. Sekerka, J. Cryst. Growth 61, 499 (1983). 
[5] G. J. Merchant and S. H. Davis, Acta Metall. 38, 2683 (1990). 

[6] W. J. Boettinger, D. Shechtman, R. J. Schaefer, and F. S. Biancaniello, Metall. Trans. A 15, 55 (1984). 

[7] M. Zimmermann, M. Carrard, and W. Kurz, Acta Metall. 37, 3305 (1989). 

[8] M. Gremaud, M. Carrard, and W. Kurz, Acta Metall. 39, 1431 (1991). 

[9] W. Kurz and R. Trivedi, Acta Metall. 38, 1 (1990). 
[10] J. J. Favier, J. Electrochem. Soc. 129, 2355 (1982). 
[11] W. van Saarloos and J. D. Weeks, Phys. Rev. Lett. 51, 1046 (1893). 
[12] B. Caroli, C. Caroli, and B. Roulet, J. Phys. (Fr.) 44, 945 (1983). 
[13] B. Caroli, C. Caroli, and B. Roulet, Acta Metall. 34, 1867 (1986). 

[14] Y. A. Buevich and V. V. Mansurov, Inzhenerno-Fizicheskii Zhurnal 47, 773 (1984), (in Russian). 
[15] Y. A. Buevich and V. V. Mansurov, Doklady AN SSSR 319, 862 (1991), (in Russian). 
[16] M. B. Geilikman and D. E. Femkin, Pisma v ZhTF 36, 238 (1982), (in Russian). 
[17] M. B. Geilikman and D. E. Femkin, J. Cryst. Growth 67, 607 (1984). 
[18] D. R. Jenkins, J. Cryst. Growth 102, 481 (1990). 

[19] V. S. Yuferev and E. N. Kolesnikova, Pisma v ZhTF 16, 76 (1991), (in Russian). 
[20] A. P. Guskov, Teplofizika Vysokikh Temperatur 29, 275 (1991), (in Russian). 

[21] W. J. Boettinger and J. H. Perepezko, in Rapid Solidified Crystalline Alloys, Proc. TMS-AIME Northeast Regional Meeting, 

edited by S. Das, B. H. Kear, and C. M. Adam (Morriston, NJ, 1985), pp. 21-58. 
[22] R. F. Wood, Appl. Phys. Lett. 37, 302 (1980). 

[23] K. A. Jackson, G. H. Gilmer, and H. J. Leamy, in Laser and Electron Beam Processing of Materials, edited by C. W. 

White and P. S. Peercy (Academic Press, New York, 1980), Chap. Solute Trapping, pp. 104-116. 
[24] M. J. Aziz, J. Appl. Phys. 53, 1158 (1982). 

[25] W. J. Boettinger and S. R. Coriell, in Rapid Solidification Materials and Technologies, edited by P. R. Sahm, H. Jones, 

and C. M. Adam (Nijhoff, Dordrecht, 1986), Chap. Microstructure Formation in Rapidly Solidified Alloys, pp. 81-108. 
[26] A. Novick-Cohen and G. I. Sivashinsky, Physica D 20, 237 (1986). 
[27] A. Novick-Cohen, Physica D 26, 403 (1987). 

[28] D. A. Huntley and S. H. Davis, Acta metallurgica et materialia, to be published . 
[29] A. Karma and A. Sarkissian, Phys. Rev. E 47, 513 (1993). 



15 



[30] F. X. Kelly and L. H. Ungar, Phys. Rev. B 34, 1746 (1986). 

[31] V. V. Voronkov, in Crystal Growth (Yerevan University, Yerevan, 1975), Vol. 11, Chap. Coefficient of Impurity Trapping 

from Melt as Function of the Orientation Angle of Stepped Surface and the Crystal Growth, pp. 357-367, (in Russian). 
[32] V. V. Voronkov, V. P. Grishin, and Y. M. Shashkov, Neorgan. Mater. 3, 2139 (1967), (in Russian). 

[33] A. A. Chernov, in Modern Crystallography (Nauka, Moscow, 1980), Vol. 3, Chap. Crystallisation Processes, p. 7, (in 
Russian) . 

[34] A. F. Witt and H. C. Gatos, J. Electrochem. Soc. 113, 808 (1966). 

[35] G. Mtiller, in Convection and Inhomogeneities in Crystal Growth from the Melt, 12. Crystals: Growth, Properties, and 

Applications, edited by H. C. Freyhardt (Springer- Verlag, Berlin Heidelberg, 1988). 
[36] R. J. Braun and S. H. Davis, J. Cryst. Growth 112, 670 (1991). 
[37] R. J. Braun and S. H. Davis, Acta Metall. 40, 2617 (1992). 
[38] R. J. Braun, G. J. Merchant, and S. H. Davis, Phys. Rev. B 45, 7002 (1992). 
[39] D. A. Huntley and S. H. Davis, J. Cryst. Growth 132, 141 (1993). 
[40] F. J. Bourland and R. Haberman, SIAM J. Appl. Math. 48, 737 (1988). 

[41] G. J. Merchant, R. J. Braun, K. Brattkus, and S. H. Davis, SIAM J. Appl. Math. 52, 1279 (1992). 
[42] K. Brattkus and D. I. Meiron, SIAM J. Appl. Math. 52, 1303 (1992). 
[43] S. M. Baer and T. Erneux, SIAM J. Appl. Math. 46, 721 (1986). 
[44] S. M. Baer and T. Erneux, SIAM J. Appl. Math. 52, 1651 (1992). 

[45] V. V. Gafiychuk, I. A. Lubashevsky, and V. V. Osipov, Sov. Phys. Crystallogr. 37 (4), 447-451 (1992). 

[46] E. Ben- Jacob, N. Goldenfeld, J. Langer, and G. Schon, Phys. Rev. A 29, 330 (1984). 

[47] M. Carrard, M. Gremaud, M. Zimmermann, and W. Kurz, Acta Metall. 40, 983 (1992). 

[48] P. Pelce, D. Pochwerger, and A. Karma, J. Cryst. Growth 110, 85 (1991). 

[49] A. Karma and N. Goldenfeld, Phys. Rev. B 31, 7018 (1985). 

[50] J. A. Warren and J. S. Langer, Phys. Rev. E 47, (1993). 

[51] N. N. Bogolyubov and Y. A. Mitropolskii, Asymptotic Methods in the Theory of Nonlinear Oscillations, 4 ed. (Nauka, 
Moscow, 1974), (in Russian). 

[52] B. D. Hassard, N. D. Kazarinoff, and Y. H. Wan, Theory and Applications of Hopf Bifurcation, 4 ed. (Cambridge University 
Press, Cambridge, 1981). 

[53] A. A. Andronov, A. A. Vitt, and S. E. Khaikin, Theory of Oscillations (Perganon Press, Oxford, 1966). 



FIG. 1. Geometry of the stability region in the CooV-plane depending on the surface tension. (The interface is stable at 
points below the solid curve. Figures (a), (b), and (c) correspond to the conditions d <SC di, d ~ di, and d 2> d±, respectively. 
Curves I and II are the boundaries of the marginal and oscillatory instabilities.) 

FIG. 2. Relative position of the absolute stability boundary (curve 1), the neutral curve (curve 2), and the curve of the 
transition between the supercritical and subcritical bifurcation (curve 3). 

FIG. 3. The dimensionless front velocity u (a) and the dimensionless interface position a (b) vs. time (in units D/V 2 ) for 
the supercritical bifurcation. 

FIG. 4. The interface oscillations and the nullcline for the supercritical bifurcation (the thick solid and dotted lines are the 
limit cycle obtained within the framework of the full diffusion problem and in the quasistationary approximation, respectively). 

FIG. 5. Distribution of the impurity concentration C so i in the growing crystal for the supercritical bifurcation (the concen- 
tration C so i and coordinate y are measured in units Coo and D/V). 



FIG. 6. Schemetic illustration of the phase path geometry in the wa-plane in limit (|4 . l| ) (the thick line represents the nullcline 
a = ao[u] and the dashed line separates the regions of closed and unclosed phase paths). 

FIG. 7. Characteristic geometry of the limit cycle in the ita-plane for the subcritical bifurcation (the solid line is the nullcline 
ao[w], the thick and dashed lines correspond to the limit cycle constracted analitically and predicted by the classic theory of 
relaxation oscillations, respectively). 



1G 



FIG. 8. The interface oscillations and the nullcline for the subcritical bifurcation (the thick and dotted lines are the limit 
cycle obtained within the framework of the full diffusion problem and in the quasistationary approximation, respectively). 

FIG. 9. The dimensionless front velocity u (a) and the dimensionless interface position a (b) vs. time (in units D/V 2 ) for 
the subcritical bifurcation. 

FIG. 10. Distribution of the impurity concentration C so \ in the growing crystal for the subcritical bifurcation (the concen- 
tration Csoi and coordinate y are measured in units Coo and D/V). 
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